N88- 17591 

EVALUATION OF A RESEARCH CIRCULATION 
CONTROL AIRFOIL USING NAVIER- STOKES METHODS 


George D. Shrevsbury 
Advanced Flight Sciences Department 
Lockheed-Georgia Company, Marietta, Georgia 


Abstract 

The compressible Reynolds time averaged Navier-Stokes equations were used 
to obtain solutions for flows about a two-dimensional circulation control 
airfoil. The governing equations were written in conservation form for a 
body-fitted coordinate system and solved using an Alternating Direction 
Implicit (ADI) procedure. A modified algebraic eddy viscosity model was 
used to define the turbulent characteristics of the flow, including the 
wall jet flow over the Coanda surface at the trailing edge. Numerical 
results are compared to experimental data obtained for a research 
circulation control airfoil geometry. Excellent agreement with the 
experimental results was obtained. 


Introduction 

One of the most efficient of the various methods for generating increased 
lift is the circulation control (CC) airfoil. This concept was developed 
in England* - and introduced into the United States by U.S. Navy 
researchers 3- . It has subsequently been the subject of extensive 
experimental test programs which have confirmed the high-lift capability of 
this innovative concept”. These airfoils obtain lift augmentation by 
tangentially exhausting a thin jet sheet over a rounded trailing edge with 
the jet sheet remaining attached well onto the airfoil lover surface due to 
the Coanda effect. 

Formerly, analysis methods for CC airfoils 11 "* -13 consisted of computational 
procedures which used weakly coupled viscous-inviscid procedures to define 
the complex flow fields resulting from the presence of the jet sheet 
exhausting into the trailing edge region. Particularly good results were 
obtained by using a potential flow CC airfoil solver developed by Dvorak, 
et al 11 ’ 13 coupled with a paraboliz^ Navier-Stokes wall-jet analysis 
program written by Dash, and associates 1 . 

The complex flow fields of the CC airfoil are governed by highly 
interactive flow regimes, however, and a comprehensive analysis of the flow 
field and the associated phenomena, including the effects of jet 
entrainment and Coanda surface geometry, requires analysis procedures which 
account for the strongly coupled nature of the viscous and inviscid flow 
regimes. 

Recently, Navier-Stokes methods have been used successfully to solve for 
the aerodynamics about CC airfoils 1 -1 . The purpose of this paper is to 
present the results obtained by using the method developed in reference 15 
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to correlate numerical results with performance data from an extensive 
experimental study 18 . This method solves the fully elliptic Navier-Stokes 
equations in 2-D planar coordinates. The mathematical and numerical 
formulations are discussed, and appropriate boundary conditions and grid 
generation procedures are defined. Modifications to existing eddy viscosity 
turbulence models to account for curved wall jets are discussed. 


Method 


Mathematical Formulation 

The development of the two-dimensional, unsteady, compressible Navier- 
Stokes equations used for this study are documented in Reference 19, and 
refinements to the method are reported in References 20-22. 

The compressible Reynolds time averaged Navier-Stokes equations may be 
written in vector form as follows: 
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The components of the viscous stress tensor are given by 
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where p is the pressure, p is the density, y is the bulk viscosity, e is 
the specific internal energy, and X is taken as -2/3 y , according to 
Stoke's hypothesis. The viscosity, y, is defined by Sutherland's lav. The 
equation of state: 

p - ptT 


is required for closure of the system of equations. 

In the above equations, all distances are normalized with respect to the 
airfoil chord, the velocities are normalized with respect to the free 
stream velocity, V®, the density is normalized with respect to the free 
stream density, and the specific internal energy is normalized with respect 
to V®^. Re and Pr are the Reynolds, number and Prandtl number, 
respectively. 


Numerical Formulation 


The numerical procedure used to solve the system of governing equations is 
a modified form of the Briley-McDonald^ Alternating Direction Implicit 
(ADI) procedure, which is based on the Douglas-Gunn^ method. It is also 
closely related to the Warming-Beam^ algorithm. Variable time steps and 
numerical dissipation have been incorporated to accelerate the convergence 
for steady state flow problems. 

The method can be outlined as follows: The governing equations are 
parabolic with respect to time. Assuming the flow field is known at a time 
level t n , the numerical procedure is used to advance the solution to a new 
time level t n+ i using a fairly large time step. If a steady state solution 
is desired, the procedure at each cell is advanced at a different time step 
based on the local cell Reynolds number. The metric terms t x , £y, etc., 
are evaluated numerically at an intermediate time level t n+ i/2- The mixed 
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derivatives that arise from terms such as (J»x^ T xx) > etc., are lagged one 
time step. The flow quantities p, u, v, and e at the new time level are 
written in terms of their values at the known time level and the 
incremental quantities; i.e., 

p n+l - p n ♦ Ap a 


The non-linear terms involved are linearized by using a Taylor expansion 
about the solution at the known time level t n . Performing these operations 
and taking all the known quantities to the right hand side, one obtains a 
linear system of equations for the incremental quantities at each grid 
point in the computational plane, excluding the boundaries. The difference 
equations may be written in matrix form as: 


(A) Uq)%3^ [i»] lCl Uq,n 
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The Douglas-Gunn procedure for generating an ADI scheme is used to solve the 
above system of equations by approximate factorization of equation (6) into 
two equations, where each involves only a one-dimensional operator: 

(A) (Aql* ♦ Jf 1#) (A q } * * (R)n (7) 
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where 
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Equations (7) and (8) are discretized using second order accurate difference 
formulas for the spatial derivatives. This technique results in a matrix 
system with a block tridiagonal structure which may be solved efficiently by 
using standard block elimination procedures. The boundary conditions for 
the unknown vector {Ag} are evaluated explicitly. Once {Ag} n is obtained, 
the flow field variables at the new time level are explicitly known. 

Artificial Dissipation Terms 

To suppress the high frequency components that appear in regions containing 
severe pressure gradients, i.e., the neighborhood of shock waves or 
stagnation points, artificial dissipation terms have been added in 
conservative form. In the present application, a blend of second and fourth 
order terms with coefficients which depend on the magnitude of the local 
pressure gradient have been added explicitly for each dependent variable in 
the manner suggested by Jameson^ , et al., and second order dissipation 
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terms have been added implicitly for each of the dependent variables. The 
coefficients of the implicit terms were added in the manner suggested by 
Steger . Extensive numerical experiments have shown that the blending of 
the dissipation terms provided better shock wave prediction with controlled 
overshoot pressure distribution. 


Turbulence Model 


28 

An algebraic eddy viscosity model developed by Baldwin and Lomax was used 
to define the turbulence transport everywhere except in the wall jet free 
shear layer. This model permits the calculation of the turbulence 
characteristics of the boundary layer by defining a two layer system. The 
viscosity in the inner layer is given by simple mixing length theory, where 
the length scale is proportional to the distance from the wall multiplied by 
the van Driest damping term, and the velocity scale is proportional to the 
length multiplied by the absolute value of the vorticity 


v - pi 2 |u| 
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where 
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In the outer layer, the velocity and length scales are constant and the 
turbulent viscosity is calculated from: 
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and yMAX is defined as the y at which FmaX occurs. F K LEB and UoiFF are 
defined by 
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The constants used in these equations are defined in Reference 28. 

The division between the inner and outer layers is taken as that point at 
which 
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The turbulence characteristics of the curved vail jet on the Coanda surface 
require special treatment, since the extra rates of strain produced by the 
curvature can exert an influence on the turbulence structure by augmenting 
or suppressing radial velocity fluctuations. In a curved vail jet, such as 
that shovn in Figure 1, a balance of centrifugal and pressure forces on a 
fluid element reveals that increases in velocity with distance from the 
center of streamline curvature generate stable flows, while flows in which 
the velocity decreases from the center of curvature are des tabilized z . In 
turbulent flow, these stabilities and instabilities lead to an increase or 
decrease in turbulent transport. This influence can result in viscosities 
which are an order of magnitude greater than those obtained in planar 
flows' 5 . Accordingly, turbulence models using standard eddy viscosity 
relations will require significant empirical modifications to reproduce the 
characteristics of curved shear layer flows. For this study, the mixing 
length was multiplied by an empirical curvature correction 

f m 1 - a S 

where a is an empirical constant whose value depends on the particular 
flow considered. A review of the literature suggests that most 
researchers place the constant in the range 6<a<14 for wall bounded flows. 
For this study, however, a value of 25 produced results more nearly in 
agreement with experimental data. S is a dimensionless parameter which is 
representative of the ratio of the extra rate of strain produced by the 
curvature to the inherent shear strain 
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where U is the velocity in the streamvise direction, n is the normal 
direction, and r is the local radius of curvature of the streamline 
considered. In areas where the curvature is small to moderate, the 
correction to the eddy viscosity is negligible. 

The location of the wake was approximated by determining the point, nearest 
the trailing edge, at which the U component of the cont ravariant velocity 
at the second grid line changed direction. The wake was then arbitrarily 
defined to exist in the region contained in the four grid points on each 
side of that location. The calculation of the eddy viscosity began at the 
lover edge of the wake, and proceeded clockwise to the upper edge of the 
wake. Values of the turbulent viscosity for the wake were then 
interpolated linearly from the values at the wake edges. 

Grid Generation and Boundary Conditions 

A body-fitted coordinate system is desired for numerical analysis 
procedures since boundary surfaces in the physical plane are mapped onto 
rectangular surfaces in the transformed plane, and the boundary conditions 
in the transformed plane may be treated more accurately. Computer methods 
developed by Thomas' 31 were employed to generate a suitable body-fitted, 
curvilinear grid system. This procedure uses a Poisson solver to define 
two-dimensional grids about airfoils and other shapes. 
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The unique trailing edge geometry and flow characteristics of CC airfoils 
makes the use of conventional C-grids difficult, since it is impossible to 
locate the cut line so that it corresponds to the physical location of the 
wake. Consequently, it was decided to use an 0-grid topology for this 
analysis. This choice represents a compromise between suitable resolution 
on the Coanda surface and adequate definition in the near-vake region. 

The grid spacing in the direction normal to the airfoil surface was 
sufficiently dense to permit satisfactory resolution of the boundary layer. 
Sixty-one grid lines were used in this direction, and approximately twenty 
of these were submerged in the boundary layer. The grid spacing in the 
normal direction varied from 0.00007 chords at the wall to 0.60 chords at 
the outer boundary. The outer boundary was defined as circular, and was 
fourteen chords in diameter. One-hundred and fifty-one points were used in 
the wrap-around direction. Grid points were clustered to permit 
satisfactory resolution at critical locations, such as the leading edge and 
blowing slot exit planes. One of the computational grids used for this 
study is shown in Figure 2. 

Boundary conditions for the computational plane consisted of specifying the 
flow conditions along the airfoil surface, including the blowing slot exit 
plane, the 0-grid cut line, and the outer boundary. On the airfoil 
surface, an adiabatic wall condition, de/dX] = 0, was imposed and 
extrapolated values of density were specified. A no-slip condition (u = v 
= 0) was used to define the velocities. At the slot blowing exit, 
specified values of total pressure and total temperature were used with an 
extrapolated value of pressure to define the boundary characteris t ics . 
Along the grid cut line, boundary conditions were applied explicitly as 
the average of the extrapolated values from each direction. 

At the outer boundary, conditions were applied according to the rule that 
flow variables should be extrapolated along characteristics leaving the 
cell and specified along characteristics entering the cell. Accordingly, 
for subsonic conditions where the boundary is experiencing inflow, values 
of the velocity and pressure are specified, while the energy is 
extrapolated from the interior. For outflow conditions, the pressure is 
specified, while values of velocity and energy are extrapolated from the 
interior. Numerical disturbances generated by the body may be reflected 
back into the computational plane, creating an adverse influence on the 
convergence characteristics of the solution. To eliminate the reflection 
of unwanted propagat ions , the pressure is specified accordingly non- 
reflecting boundary criteria prescribed by Rudy and Strikwerda , which 
have been implemented at the outflow boundaries. 

In all cases where extrapolated values were specified at the boundaries, a 
two-point extrapolation of the form 

4 1 

q l " I q 2 * 3 q 3 


was used. 


Results 
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Research CC Airfoil 


1 ft 

Novak and Cornelius conducted wind tunnel tests on a 15.6 per cent thick 
CC airfoil section which had been specifically designed to provide data for 
Navier-Stokes code validation. The bloving slot height-to-radius ratio was 
0.1, and the overall chord length was 15 inches. This model was designed 
with a cylindrical Coanda surface with a radius-to-chord ratio of 0.067. 
While this ratio is relatively high and is certainly not representative of 
practical flight systems, it does provide a physically large slot height, 
which improves the quality of the measurements. Data were acquired in the 
Lockheed-Georgia Low Turbulence Wind Tunnel at a free-stream Mach number of 
0.0853 and a Reynolds number of 780,000. The model angle of attack was 
zero degrees. The data consisted of airfoil surface pressure measurements 
and extensive flow field surveys using a laser velocimeter (LV). The 
general profile of the section can be visualized from the grid shown in 
Figure 2. 

Once a suitable computaional grid had been constructed, numerical studies 
were conducted at a Mach number and Reynolds number corresponding to the 
experimental tests. The angle of attack was varied numerically until a 
lift coefficient based on integrated pressures was obtained which 
corresponded to the experimental zero incidence case. For the jet total 
pressure ratio invest igated , the numerical angle of attack was -2 degrees, 
and the corresponding lift coefficient was 4.55. All computed and 
experimental data subsequently presented are for that lift coefficient. 
The jet total pressure ratio was 1.10, and the ratio of jet total 
temperature to free-stream total temperature was .964. 

The numerical results were obtained by executing the code on the 
Lockheed/ASG Cray X-MP/24 computer. Approximately 1000 iterations were 
required to obtain a converged, steady-state solution. This formulation 
of the Navier-Stokes equations requires approximately 2.5 X 10 CPU 
seconds/grid-point/ time-step of Cray execution time. 

Computed streamlines for the research airfoil are shown in Figure 3. This 
figure clearly demonstrates the characteristic flows for CC airfoils, 
including the large, induced circulation which produces a strong downwash 
in the wake region and upwash at the leading edge. The jet entrainment 
effects on the upper surface and the Coanda turning of the jet can also be 
seen. 

Computed velocity vectors and streamlines in the trailing edge region are 
shown for the same case in Figure 4. Details of wall jet development, as 
well as the interaction of the upper and lower surface flows can be clearly 
visualized. These results demonstrate the attached, well-behaved nature of 
the flow, which is characteristic of CC airfoils, even near stall 5 . 

Comparisons of computed and experimental velocity profiles at the jet exit 
plane are shown in Figure 5. The profiles are shown as ratios of the local 
velocity to the free-stream velocity versus the radial distance above the 
Coanda surface. The upper edge of the exhaust jet is located at a Y of .1 
inches. The data in Figure 5(a) include the boundary layer which has been 
established on the airfoil upper surface as well as the jet slot exhaust 
flow. The comparison between the computed and experimental boundary layers 
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is very good, which suggests that the jet entrainment effects propagated 
upstream are being properly modeled. 

A detailed comparison of the computed and experimental velocity profiles 
for the jet exhaust is shown in Figure 5(b). In order to produce 
corresponding values of jet exhaust velocities, the numerical data were 
run at a jet total pressure ratio of 1.10, compared to a value of 1.12 
measured experimentally. The discrepancy may be due to total pressure 
losses experienced in the duct between the plenum measurement location and 
the jet exit plane. The numerical total pressure ratio used however, 
provides an excellent reproduction of the experimental velocity profile 
except at the upper edge of the jet, where an established jet boundary 
layer already existed. The consequences of failing to properly model this 
characteristic of the jet exhaust profile are not known. 

The jet exit total pressure was assumed to be constant across the jet slot 
height except at the walls, where a no-slip ( u = v = 0) condition was 
enforced. The small velocity deficit occurring near the bottom of the 
numerical profile is believed to be the result of errors introduced by 
poor grid characteristics in that region. The comparison between the 
computed and experimental velocity gradients produced by the strong radial 
pressure differentials is very good. The numerical velocities obtained 
from this total pressure ratio resulted in a jet momentum coefficient of 
approximaely 0.275. 

It is interesting to note that the differences in experimental boundary 
layers observed at the upper and lower edges of the jet are consistent with 
the previously discussed postulate that positive radial velocity gradients 
in regions of curvature are stabilizing, while negative velocity gradients 
in that direction are destabilizing. 

Experimental and computed velocity profiles for the wall jet flow at a 
circumferencial location of 90 degrees, measured clockwise from the jet 
exit location, are- shown in Figure 6. The conditions at which these 
results were obtained are the same as for the previous figure. While the 
general magnitudes and wall jet thicknesses agree reasonably well, the 
differences in profile characteristics near the wall are significant. The 
extremely stable characteris tics exhibited by the experimental flow 
adjacent to the wall are not reproduced adequately by the numerical data. 
Since the experimental conditions at the beginning of the Coanda region are 
closely approximated, it is concluded that an empirically corrected eddy 
viscosity turbulence formulation is not sufficient to properly model the 
turbulent characterist ics of curved wall jets, particularly the strong 
stabilization that occurs near the boundary. As a consequence of this 
discrepancy, the numerical wall jet dissipates energy too rapidly and 
experiences premature detachment. The angle of attack correction of -2 
degrees was underpredicted therefore, and the actual equivalent 
experimental angle of attack was somewhat more negative. 

i 

Experimental and computed CC airfoil pressure distributions are shown in 
Figure 7. The agreement between experimental and computed pressure 
distributions is very good. The strong suction peaks produced by the 
super-circulation at the leading edge and the jet sheet turning on the 
Coanda surface are very accurately predicted. The discrepancy between the 
computed and experimental data on the lower surface, near the trailing 
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edge, is probably the result of differences in the locations of the jet 
sheet detachment points, as observed experimentally and predicted 
numerically. 


Conclusions 

A computational procedure has been developed which permits the calculation 
of the performance characteristics of circulation control airfoils over a 
broad range of free-stream conditions. The fully elliptic, Reynolds time 
averaged Navier-Stokes equations were solved numerically, using an 
Alternating Direction Implicit (ADI) algorithm. The computed results 
compared well with experiments conducted on a research CC airfoil which had 
been specifically designed to provide data for Navier-Stokes code 
validation, including force data and detailed flow measurements taken in 
the trailing edge region. A specially modified algebraic eddy viscosity 
model was used to predict the behavior of the wall jet, and although the 
overall behavior of the curved wall jet was sufficiently approximated, 
important turbulent characteristics crucial to the prediction of the jet 
sheet detachment point were not adequately predicted. Extensions to the 
present work will include the incorporation of advanced turbulence models 
to provide improved analysis of the wall jet characteristics. 
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Figure 2. Computational Grid 
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Figure 3. Computed Streamlines for Research CC Airfoil 
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Figure 4. Computed Ve] 












x/c 

Figure 7. Computed and Experimental Pressure Coefficient 
Distributions 
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